//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : RrbSolver.h
//----------------------------------------------------------------------------

#ifndef RRBSOLVER_H
#define RRBSOLVER_H

#include "Solver.h"
#include "Array2D.h"

#define RRBSOLVER_BASE Solver

//----------------------------------------------------------------------------
// class RrbSolver
//----------------------------------------------------------------------------
/// recursive red black solver by Auke van der Ploeg
//----------------------------------------------------------------------------
class RrbSolver : public RRBSOLVER_BASE
{
    typedef RRBSOLVER_BASE base;

public:

    RrbSolver(Array2D<REAL> *pcc,
              Array2D<REAL> *pcs,
              Array2D<REAL> *pcw,
              int Nx1,
              int Nx2,
              REAL tolerance);

    virtual ~RrbSolver();

    virtual void solve(Array2D<REAL> *pB,
                       Array2D<REAL> *pX);

private:

    /// @cond nevercall
    RrbSolver(const RrbSolver& m);
    RrbSolver & operator =(const RrbSolver& m);
    /// @endcond

    void preconditioner();

    void solver();

    void matvec(const Array2D<REAL> &X, Array2D<REAL> *Az);

    Array2D<REAL> *m_pcc; /// pointer to MatC
    Array2D<REAL> *m_pcs; /// pointer to MatS
    Array2D<REAL> *m_pcw; /// pointer to MatW
    Array2D<REAL> m_cne; /// K1
    Array2D<REAL> m_cnw; /// K2
    Array2D<REAL> m_cse; /// Az
    Array2D<REAL> m_csw; /// K4
    Array2D<REAL> m_ccc; /// copyMatC
    Array2D<REAL> m_css; /// copyMatS
    Array2D<REAL> m_cww; /// copyMatW
    Array2D<REAL> m_cee; /// copyMatE
    Array2D<REAL> m_cnn; /// copyMatN

    int m_st;
    int m_tw;
    int m_niv;
    REAL m_rhobew;
    REAL m_bet;
    REAL m_alf;
    bool m_done;
    bool m_done1;
    bool m_done2;
};

inline Solver* NewRrbSolver(Array2D<REAL> *pmatC, ///< central matrix
                            Array2D<REAL> *pmatS, ///< south matrix
                            Array2D<REAL> *pmatW, ///< west matrix
                            int N, ///< width
                            int M, ///< height
                            REAL tolerance) ///< solver tolerance
{
    return new RrbSolver(pmatC, pmatS, pmatW, N, M, tolerance);
}

#endif // RRBSOLVER_H
